Moved from Explore Processed UVP so that there is less going on per notebook.

Rerun that notebook so we have the data we need.

library(tidyverse)
library(cowplot)
library(plotly)

Load all data

unbinned_DepthSummary <- read_csv("dataOut/unbinned_DepthSummary.csv")
Parsed with column specification:
cols(
  .default = col_double(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.
unbinned_EachSize <- read_csv("dataOut/unbinned_EachSize.csv")
Parsed with column specification:
cols(
  profile = col_character(),
  time = col_datetime(format = ""),
  depth = col_double(),
  psd_gam = col_double(),
  vol = col_double(),
  sizeclass = col_character(),
  lb = col_double(),
  ub = col_double(),
  binsize = col_double(),
  TotalParticles = col_double(),
  nparticles = col_double(),
  n_nparticles = col_double(),
  biovolume = col_double(),
  speed = col_double(),
  flux = col_double(),
  flux_fit = col_double(),
  GamPredictTP = col_double()
)
binned_DepthSummary <- read_csv("dataOut/binned_DepthSummary.csv")
Parsed with column specification:
cols(
  .default = col_double(),
  profile = col_character(),
  time = col_datetime(format = "")
)
See spec(...) for full column specifications.
binned_EachSize <- read_csv("dataOut/binned_EachSize.csv")
Parsed with column specification:
cols(
  profile = col_character(),
  time = col_datetime(format = ""),
  depth = col_double(),
  psd_gam = col_double(),
  lb = col_double(),
  ub = col_double(),
  binsize = col_double(),
  vol = col_double(),
  TotalParticles = col_double(),
  nparticles = col_double(),
  n_nparticles = col_double(),
  biovolume = col_double(),
  speed = col_double(),
  flux = col_double(),
  flux_fit = col_double(),
  GamPredictTP = col_double()
)
#ksource("Explore_Processed_UVP.Rmd")

Confedence interevals and PSD

Do slopes fall outside of predicted confidence intervals assuming a poisson distribution?

Unbinned

Just look at station 43, which is the deepest cast and stops the least frequently.

unbinned043 <- unbinned_DepthSummary %>% filter(profile == "stn_043") %>% arrange(depth)
plot_psd_2 <- unbinned_DepthSummary %>% ggplot(aes(x = psd, y = depth, colour = profile)) + geom_point(alpha = 0.5, size = 2, shape = 21, fill = "gray") + scale_y_reverse()  + guides(colour = FALSE) + labs(x = "Particle Size Distribution Slope", y = "Depth (m)") + theme_cowplot() + 
  geom_path(aes(x = qpt05, y = depth), color = "blue", data = unbinned043) +
  geom_path(aes(x = qpt95, y = depth), color = "blue", data = unbinned043) +
  geom_path(aes(x = psd_gam, y = depth), color = "black", data = unbinned043)
ggplotly(plot_psd_2)

There’s a little noise in the quantiles, which I think is because the volumes change from one bin to the nest.

Binned

As above, but binned. Only looking at station 43 here.

binned043 <- binned_DepthSummary %>% filter(profile == "stn_043") %>% arrange(depth)
plot_psd_2b <- binned_DepthSummary %>%
  ggplot(aes(x = psd, y = depth, colour = profile)) + geom_point(alpha = 0.5, size = 2, shape = 21, fill = "black") +
  scale_y_reverse(breaks = seq(from = 0, to = 2500, by = 500))  + guides(colour = FALSE) + labs(x = "Particle Size Distribution Slope", y = "Depth (m)") + theme_cowplot() + 
  geom_path(aes(x = qpt05, y = depth), data = binned043, color = "blue")  + geom_path(aes(x = qpt95, y = depth), data = binned043, color = "blue")
#ggplotly(plot_psd_2b)
plot_psd_2b

Even with binning there are many points that fall outside of the gam. So binning at this level doesn’t really slay the variance.

Upper size cutoff

Here we ask whether there is an upper size cutoff to particle size, that we can detect.

For one/each depth, identify the non-zero points, fit a gam to those. Then see if we would predict zero particles in the next smallist bin.

If we would expect to see particles, but we don’t something is removing the bigger than that particles.

binned_DepthSummary %>% filter(profile == "stn_043", depth == 275)
test_depth <- binned_EachSize %>% filter(profile == "stn_043", depth == 275)
test_depth
fit_model = function(df) glm(TotalParticles ~ log(lb), offset = log(binsize * vol), data = df, family = "poisson")
test_depth_sawParticles <- test_depth %>% filter(TotalParticles >= 1)
test_depth_firstZero <- test_depth %>% filter(TotalParticles == 0) %>% filter(lb == max(lb))
test_mod <- fit_model(test_depth_sawParticles)
pred_firstZero <- predict(test_mod, test_depth_firstZero, type = "response")
lower_pt05 <- qpois(0.05, lambda = pred_firstZero)
pred_firstZero
           1 
8.102217e-05 
lower_pt05
[1] 0

So in this particular case, the first zero is totally expected. Are there exceptions?

expected_firstZero <- function(df){
  loc_depth_sawParticles <- df %>% filter(TotalParticles >= 1)
  loc_depth_firstZero <- df %>% filter(TotalParticles == 0) %>% filter(lb == max(lb))
  loc_mod <- fit_model(loc_depth_sawParticles)
  loc_pred_firstZero <- predict(test_mod, test_depth_firstZero, type = "response")
  loc_lower_pt05 <- qpois(0.05, lambda = pred_firstZero)
  loc_lower_pt05
}

expected_firstZero(test_depth)
[1] 0
look_nz <- binned_EachSize %>% group_by(profile, time, depth) %>% nest() %>% 
  mutate(fz = map_dbl(data, expected_firstZero))
look_nz
ggplot(look_nz, aes(y = depth, x = fz)) + geom_point() + scale_y_reverse()

Ok. So, um. There don’t appear to be cases where we would expect to see a non zero where we actually see a zero. This says to me that UVP is not a suitable way to define upper bounds of particle sizes.

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKCk1vdmVkIGZyb20gRXhwbG9yZSBQcm9jZXNzZWQgVVZQIHNvIHRoYXQgdGhlcmUgaXMgbGVzcyBnb2luZyBvbiBwZXIgbm90ZWJvb2suCgpSZXJ1biB0aGF0IG5vdGVib29rIHNvIHdlIGhhdmUgdGhlIGRhdGEgd2UgbmVlZC4KCmBgYHtyLCBtZXNzYWdlID0gRkFMU0V9CmxpYnJhcnkodGlkeXZlcnNlKQpsaWJyYXJ5KGNvd3Bsb3QpCmxpYnJhcnkocGxvdGx5KQpgYGAKCiMgTG9hZCBhbGwgZGF0YQpgYGB7cn0KdW5iaW5uZWRfRGVwdGhTdW1tYXJ5IDwtIHJlYWRfY3N2KCJkYXRhT3V0L3VuYmlubmVkX0RlcHRoU3VtbWFyeS5jc3YiKQp1bmJpbm5lZF9FYWNoU2l6ZSA8LSByZWFkX2NzdigiZGF0YU91dC91bmJpbm5lZF9FYWNoU2l6ZS5jc3YiKQpiaW5uZWRfRGVwdGhTdW1tYXJ5IDwtIHJlYWRfY3N2KCJkYXRhT3V0L2Jpbm5lZF9EZXB0aFN1bW1hcnkuY3N2IikKYmlubmVkX0VhY2hTaXplIDwtIHJlYWRfY3N2KCJkYXRhT3V0L2Jpbm5lZF9FYWNoU2l6ZS5jc3YiKQpgYGAKCgpgYGB7ciwgbWVzc2FnZSA9IEZBTFNFfQoja3NvdXJjZSgiRXhwbG9yZV9Qcm9jZXNzZWRfVVZQLlJtZCIpCmBgYAoKCiMgQ29uZmVkZW5jZSBpbnRlcmV2YWxzIGFuZCBQU0QKRG8gc2xvcGVzIGZhbGwgb3V0c2lkZSBvZiBwcmVkaWN0ZWQgY29uZmlkZW5jZSBpbnRlcnZhbHMgYXNzdW1pbmcgYSBwb2lzc29uIGRpc3RyaWJ1dGlvbj8KCiMjIFVuYmlubmVkCgpKdXN0IGxvb2sgYXQgc3RhdGlvbiA0Mywgd2hpY2ggaXMgdGhlIGRlZXBlc3QgY2FzdCBhbmQgc3RvcHMgdGhlIGxlYXN0IGZyZXF1ZW50bHkuCmBgYHtyfQp1bmJpbm5lZDA0MyA8LSB1bmJpbm5lZF9EZXB0aFN1bW1hcnkgJT4lIGZpbHRlcihwcm9maWxlID09ICJzdG5fMDQzIikgJT4lIGFycmFuZ2UoZGVwdGgpCmBgYAoKCmBgYHtyfQpwbG90X3BzZF8yIDwtIHVuYmlubmVkX0RlcHRoU3VtbWFyeSAlPiUgZ2dwbG90KGFlcyh4ID0gcHNkLCB5ID0gZGVwdGgsIGNvbG91ciA9IHByb2ZpbGUpKSArIGdlb21fcG9pbnQoYWxwaGEgPSAwLjUsIHNpemUgPSAyLCBzaGFwZSA9IDIxLCBmaWxsID0gImdyYXkiKSArIHNjYWxlX3lfcmV2ZXJzZSgpICArIGd1aWRlcyhjb2xvdXIgPSBGQUxTRSkgKyBsYWJzKHggPSAiUGFydGljbGUgU2l6ZSBEaXN0cmlidXRpb24gU2xvcGUiLCB5ID0gIkRlcHRoIChtKSIpICsgdGhlbWVfY293cGxvdCgpICsgCiAgZ2VvbV9wYXRoKGFlcyh4ID0gcXB0MDUsIHkgPSBkZXB0aCksIGNvbG9yID0gImJsdWUiLCBkYXRhID0gdW5iaW5uZWQwNDMpICsKICBnZW9tX3BhdGgoYWVzKHggPSBxcHQ5NSwgeSA9IGRlcHRoKSwgY29sb3IgPSAiYmx1ZSIsIGRhdGEgPSB1bmJpbm5lZDA0MykgKwogIGdlb21fcGF0aChhZXMoeCA9IHBzZF9nYW0sIHkgPSBkZXB0aCksIGNvbG9yID0gImJsYWNrIiwgZGF0YSA9IHVuYmlubmVkMDQzKQpnZ3Bsb3RseShwbG90X3BzZF8yKQpgYGAKVGhlcmUncyBhIGxpdHRsZSBub2lzZSBpbiB0aGUgcXVhbnRpbGVzLCB3aGljaCBJIHRoaW5rIGlzIGJlY2F1c2UgdGhlIHZvbHVtZXMgY2hhbmdlIGZyb20gb25lIGJpbiB0byB0aGUgbmVzdC4KCiMjIEJpbm5lZAoKQXMgYWJvdmUsIGJ1dCBiaW5uZWQuIE9ubHkgbG9va2luZyBhdCBzdGF0aW9uIDQzIGhlcmUuCmBgYHtyfQpiaW5uZWQwNDMgPC0gYmlubmVkX0RlcHRoU3VtbWFyeSAlPiUgZmlsdGVyKHByb2ZpbGUgPT0gInN0bl8wNDMiKSAlPiUgYXJyYW5nZShkZXB0aCkKcGxvdF9wc2RfMmIgPC0gYmlubmVkX0RlcHRoU3VtbWFyeSAlPiUKICBnZ3Bsb3QoYWVzKHggPSBwc2QsIHkgPSBkZXB0aCwgY29sb3VyID0gcHJvZmlsZSkpICsgZ2VvbV9wb2ludChhbHBoYSA9IDAuNSwgc2l6ZSA9IDIsIHNoYXBlID0gMjEsIGZpbGwgPSAiYmxhY2siKSArCiAgc2NhbGVfeV9yZXZlcnNlKGJyZWFrcyA9IHNlcShmcm9tID0gMCwgdG8gPSAyNTAwLCBieSA9IDUwMCkpICArIGd1aWRlcyhjb2xvdXIgPSBGQUxTRSkgKyBsYWJzKHggPSAiUGFydGljbGUgU2l6ZSBEaXN0cmlidXRpb24gU2xvcGUiLCB5ID0gIkRlcHRoIChtKSIpICsgdGhlbWVfY293cGxvdCgpICsgCiAgZ2VvbV9wYXRoKGFlcyh4ID0gcXB0MDUsIHkgPSBkZXB0aCksIGRhdGEgPSBiaW5uZWQwNDMsIGNvbG9yID0gImJsdWUiKSAgKyBnZW9tX3BhdGgoYWVzKHggPSBxcHQ5NSwgeSA9IGRlcHRoKSwgZGF0YSA9IGJpbm5lZDA0MywgY29sb3IgPSAiYmx1ZSIpCiNnZ3Bsb3RseShwbG90X3BzZF8yYikKcGxvdF9wc2RfMmIKYGBgCgoKRXZlbiB3aXRoIGJpbm5pbmcgdGhlcmUgYXJlIG1hbnkgcG9pbnRzIHRoYXQgZmFsbCBvdXRzaWRlIG9mIHRoZSBnYW0uIFNvIGJpbm5pbmcgYXQgdGhpcyBsZXZlbCBkb2Vzbid0IHJlYWxseSBzbGF5IHRoZSB2YXJpYW5jZS4KCiMgVXBwZXIgc2l6ZSBjdXRvZmYKCkhlcmUgd2UgYXNrIHdoZXRoZXIgdGhlcmUgaXMgYW4gdXBwZXIgc2l6ZSBjdXRvZmYgdG8gcGFydGljbGUgc2l6ZSwgdGhhdCB3ZSBjYW4gZGV0ZWN0LgoKRm9yIG9uZS9lYWNoIGRlcHRoLCBpZGVudGlmeSB0aGUgbm9uLXplcm8gcG9pbnRzLCBmaXQgYSBnYW0gdG8gdGhvc2UuClRoZW4gc2VlIGlmIHdlIHdvdWxkIHByZWRpY3QgemVybyBwYXJ0aWNsZXMgaW4gdGhlIG5leHQgc21hbGxpc3QgYmluLgoKSWYgd2Ugd291bGQgZXhwZWN0IHRvIHNlZSBwYXJ0aWNsZXMsIGJ1dCB3ZSBkb24ndCBzb21ldGhpbmcgaXMgcmVtb3ZpbmcgdGhlIGJpZ2dlciB0aGFuIHRoYXQgcGFydGljbGVzLgoKYGBge3J9CmJpbm5lZF9EZXB0aFN1bW1hcnkgJT4lIGZpbHRlcihwcm9maWxlID09ICJzdG5fMDQzIiwgZGVwdGggPT0gMjc1KQp0ZXN0X2RlcHRoIDwtIGJpbm5lZF9FYWNoU2l6ZSAlPiUgZmlsdGVyKHByb2ZpbGUgPT0gInN0bl8wNDMiLCBkZXB0aCA9PSAyNzUpCnRlc3RfZGVwdGgKYGBgCgoKYGBge3J9CmZpdF9tb2RlbCA9IGZ1bmN0aW9uKGRmKSBnbG0oVG90YWxQYXJ0aWNsZXMgfiBsb2cobGIpLCBvZmZzZXQgPSBsb2coYmluc2l6ZSAqIHZvbCksIGRhdGEgPSBkZiwgZmFtaWx5ID0gInBvaXNzb24iKQp0ZXN0X2RlcHRoX3Nhd1BhcnRpY2xlcyA8LSB0ZXN0X2RlcHRoICU+JSBmaWx0ZXIoVG90YWxQYXJ0aWNsZXMgPj0gMSkKdGVzdF9kZXB0aF9maXJzdFplcm8gPC0gdGVzdF9kZXB0aCAlPiUgZmlsdGVyKFRvdGFsUGFydGljbGVzID09IDApICU+JSBmaWx0ZXIobGIgPT0gbWF4KGxiKSkKdGVzdF9tb2QgPC0gZml0X21vZGVsKHRlc3RfZGVwdGhfc2F3UGFydGljbGVzKQpwcmVkX2ZpcnN0WmVybyA8LSBwcmVkaWN0KHRlc3RfbW9kLCB0ZXN0X2RlcHRoX2ZpcnN0WmVybywgdHlwZSA9ICJyZXNwb25zZSIpCmxvd2VyX3B0MDUgPC0gcXBvaXMoMC4wNSwgbGFtYmRhID0gcHJlZF9maXJzdFplcm8pCnByZWRfZmlyc3RaZXJvCmxvd2VyX3B0MDUKYGBgCgpTbyBpbiB0aGlzIHBhcnRpY3VsYXIgY2FzZSwgdGhlIGZpcnN0IHplcm8gaXMgdG90YWxseSBleHBlY3RlZC4gQXJlIHRoZXJlIGV4Y2VwdGlvbnM/CgpgYGB7cn0KZXhwZWN0ZWRfZmlyc3RaZXJvIDwtIGZ1bmN0aW9uKGRmKXsKICBsb2NfZGVwdGhfc2F3UGFydGljbGVzIDwtIGRmICU+JSBmaWx0ZXIoVG90YWxQYXJ0aWNsZXMgPj0gMSkKICBsb2NfZGVwdGhfZmlyc3RaZXJvIDwtIGRmICU+JSBmaWx0ZXIoVG90YWxQYXJ0aWNsZXMgPT0gMCkgJT4lIGZpbHRlcihsYiA9PSBtYXgobGIpKQogIGxvY19tb2QgPC0gZml0X21vZGVsKGxvY19kZXB0aF9zYXdQYXJ0aWNsZXMpCiAgbG9jX3ByZWRfZmlyc3RaZXJvIDwtIHByZWRpY3QodGVzdF9tb2QsIHRlc3RfZGVwdGhfZmlyc3RaZXJvLCB0eXBlID0gInJlc3BvbnNlIikKICBsb2NfbG93ZXJfcHQwNSA8LSBxcG9pcygwLjA1LCBsYW1iZGEgPSBwcmVkX2ZpcnN0WmVybykKICBsb2NfbG93ZXJfcHQwNQp9CgpleHBlY3RlZF9maXJzdFplcm8odGVzdF9kZXB0aCkKYGBgCgoKYGBge3J9Cmxvb2tfbnogPC0gYmlubmVkX0VhY2hTaXplICU+JSBncm91cF9ieShwcm9maWxlLCB0aW1lLCBkZXB0aCkgJT4lIG5lc3QoKSAlPiUgCiAgbXV0YXRlKGZ6ID0gbWFwX2RibChkYXRhLCBleHBlY3RlZF9maXJzdFplcm8pKQpsb29rX256CmBgYAoKYGBge3J9CmdncGxvdChsb29rX256LCBhZXMoeSA9IGRlcHRoLCB4ID0gZnopKSArIGdlb21fcG9pbnQoKSArIHNjYWxlX3lfcmV2ZXJzZSgpCmBgYAoKT2suIFNvLCB1bS4gVGhlcmUgZG9uJ3QgYXBwZWFyIHRvIGJlIGNhc2VzIHdoZXJlIHdlIHdvdWxkIGV4cGVjdCB0byBzZWUgYSBub24gemVybyB3aGVyZSB3ZSBhY3R1YWxseSBzZWUgYSB6ZXJvLgpUaGlzIHNheXMgdG8gbWUgdGhhdCBVVlAgaXMgbm90IGEgc3VpdGFibGUgd2F5IHRvIGRlZmluZSB1cHBlciBib3VuZHMgb2YgcGFydGljbGUgc2l6ZXMu